#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <cuda.h>

//funcao para calcular o valor absoluto (modulo) de um numero
float modulo(float valor);
//funcao para verificar se o sistema converge pelo criterio das linhas
int convergenciaLinha(float **novaMA, int J_ORDER);
//funcao para verificar se o sistema converge pelo criterio das colunas
int convergenciaColuna(float **novaMA, int J_ORDER);

/*
//funcao para ser executada no device
//calcula os X[i] de cada iteracao
__global__ void calcx(int J_ORDER, float* dev_novaMA, float* dev_B, float* dev_X, float* dev_novoX, float* dev_maxDif, float* dev_maxnovoX){
	int j, k;
	float modulo=0;
	float somatorio = 0;
	int tid = blockIdx.x;
	
	for(j = 0; j < J_ORDER; j++){
		if(tid != j){
			k = (tid*J_ORDER) + j;
			somatorio = somatorio + (dev_novaMA[k]*dev_X[j]);
		}
	}
	//calcula o novo valor de X
	dev_novoX[tid] = dev_B[tid]-somatorio;
	//calcula o maxDif
	modulo = dev_novoX[tid]-dev_X[tid];
	if(modulo >= 0)
		dev_maxDif[tid] = modulo;
	else
		dev_maxDif[tid] = -modulo;
	//calcula o maxnovoX
	modulo =dev_novoX[tid];
	if(modulo >= 0)
		dev_maxnovoX[tid] = modulo;
	else
		dev_maxnovoX[tid] = -modulo;

}
*/

//resolve<<<J_ORDER,1>>>(J_ORDER, dev_novaMA, dev_B, dev_novoX, J_ERROR, J_ITE_MAX);
__global__ void resolve(int J_ORDER, float *dev_novaMA, float *dev_B, float *dev_X, float *dev_novoX, float J_ERROR, int J_ITE_MAX, float *dev_maxDif, float *dev_maxnovoX)
{
    float Mr = 100.0f;
    int iteracoes = 1;//a iteracao zero eh feita manualmente

    float maxDif, maxnovoX;//variaveis auxiliares para o calculo do erro Mr

	//se passou em algum teste, realiza demais iteracoes
    while(Mr > J_ERROR && iteracoes < J_ITE_MAX)
    {//iteracoes
	maxDif = 0.0f;
        maxnovoX = 0.0f;
			
	//calcula proximo X
	//calcx(J_ORDER, dev_novaMA, dev_B, dev_X, dev_novoX, dev_maxDif, dev_maxnovoX);
	int j, k;
	float modulo = 0.0f;
	float somatorio = 0.0f;
	int tid = blockIdx.x;
	
	for(j = 0; j < J_ORDER; j++)
	{
		if(tid != j)
		{
			k = (tid*J_ORDER) + j;
			somatorio = somatorio + (dev_novaMA[k]*dev_X[j]);
		}
	}
	//calcula o novo valor de X
	dev_novoX[tid] = dev_B[tid]-somatorio;
	printf("dev_novoX[tid=%d] (%f) = dev_B[tid=%d] (%f) - somatorio (%f)\n", 
	       iteracoes, tid, dev_novoX[tid], tid, dev_B[tid], somatorio);
	//calcula o maxDif
	modulo = dev_novoX[tid]-dev_X[tid];
	if(modulo >= 0.0f)
		dev_maxDif[tid] = modulo;
	else
		dev_maxDif[tid] = -modulo;
	//calcula o maxnovoX
	modulo =dev_novoX[tid];
	if(modulo >= 0.0f)
		dev_maxnovoX[tid] = modulo;
	else
		dev_maxnovoX[tid] = -modulo;
	
/*
			cudaMemcpy(X, dev_X, J_ORDER*sizeof(float), cudaMemcpyDeviceToHost);
			cudaMemcpy(novoX, dev_novoX, J_ORDER*sizeof(float), cudaMemcpyDeviceToHost);
			cudaMemcpy(vmaxDif, dev_maxDif, J_ORDER*sizeof(float), cudaMemcpyDeviceToHost);
			cudaMemcpy(vmaxnovoX, dev_maxnovoX, J_ORDER*sizeof(float), cudaMemcpyDeviceToHost);
*/
	dev_x[tid] = dev_novoX[tid];
	__syncthreads();
	int i;		
	//calcula o maxDif
	//calcula o maxnovoX
	for(i = 0; i < J_ORDER; i++){
		if(dev_maxDif[i] > maxDif)
			maxDif = dev_maxDif[i];
		if(dev_maxnovoX[i] > maxnovoX)
			maxnovoX = dev_maxnovoX[i];
		//atualiza o valor de X[]
		//dev_X[i] = dev_novoX[i];
	}
		
	__syncthreads();
        //calcula o erro
        Mr = maxDif/maxnovoX;
	printf("[it=%d, tid=%d]Mr = %f\n", iteracoes, tid, Mr);

/*
			//atualiza o valor de X no device
			cudaMemcpy(dev_X, X, J_ORDER*sizeof(float), cudaMemcpyHostToDevice);
*/
		//vai para a proxima iteracao
 		iteracoes++;
	}//fim das iteracoes

}


int main(int argc, char **argv)
{
    /*declaracao de variaveis*/
    int J_ORDER, J_ROW_TEST, J_ITE_MAX;//ordem, fila e numero max de iteracoes
    float J_ERROR;//erro permitido
    float **MA;//matriz A
    float **novaMA;//matriz A
    float *B;//vetor B
    float *X;//armazena os X, apenas da iteracao atual
    float *novoX;//armazena os novoX, apenas da iteracao atual
    float Mr = 100;//erro comecando com um valor alto para passar na primeira iteracao
    float maxDif, maxnovoX;//variaveis auxiliares para o calculo do erro Mr
	float *vmaxDif, *vmaxnovoX;//vetores auxiliares para o calculo do erro Mr
    //float somatorio;//variavel auxiliar para substituir o somatorio no calculo do novoX
    float resultado = 0;//resultado do teste final com o J_ROW_TEST
    float resultadoEsperado;
    char* arquivo;//variável com o nome do arquivo
    FILE *pFile;//variavel para arquivo
    int i, j, k;//contadores auxiliares
    int iteracoes = 1;//a iteracao zero eh feita manualmente
    //arquivo = "exemplo.txt";//nome do arquivo
	/*variaveis para executar no device do cuda*/
	float *novaMA2;
	float *dev_novaMA;
	float *dev_B;
	float *dev_X;
	float *dev_novoX;
	float *dev_maxDif;
	float *dev_maxnovoX;
	
	//verifica os argumentos passados
	if(argc != 2){
		printf("\nQuantidade errada de argumentos..\n");
		return -1;
	}
	else
		arquivo = argv[1];
	
	//abrindo arquivo para leitura
	pFile = fopen(arquivo, "r");

    //Pegando variaveis do arquivo
    fscanf(pFile, "%i", &J_ORDER);//ordem da matriz
    fscanf(pFile, "%i", &J_ROW_TEST);//fila para ser avaliada
    fscanf(pFile, "%f", &J_ERROR);//erro permitido
    fscanf(pFile, "%i", &J_ITE_MAX);//numero maximo de iteracoes


	
    //alocando espaço para os vetores e matrizes
    MA = (float **) calloc(J_ORDER, sizeof(float*));//espaco para a matriz A
    for(i=0; i<J_ORDER; i++){
        MA[i] = (float *) calloc(J_ORDER, sizeof(float));
    }
    novaMA = (float **) calloc(J_ORDER, sizeof(float*));//espaco para a matriz A*
    for(i=0; i<J_ORDER; i++){
        novaMA[i] = (float *) calloc(J_ORDER, sizeof(float));
    }
    B = (float *) calloc(J_ORDER, sizeof(float));//espaco para o vetor B
    X = (float *) calloc(J_ORDER, sizeof(float));//espaco para o vetor X
    novoX = (float *) calloc(J_ORDER, sizeof(float));//espaco para o vetor novoX
	novaMA2 = (float *) calloc(J_ORDER*J_ORDER, sizeof(float));//espaco para o vetor auxiliar novaMA2
	vmaxDif = (float *) calloc(J_ORDER, sizeof(float));//
	vmaxnovoX = (float *) calloc(J_ORDER, sizeof(float));//
	

	
    //pegando os valores da matriz A
    for(i=0;i<J_ORDER;i++){
        for(j=0;j<J_ORDER;j++){
            fscanf(pFile, "%f", &MA[i][j]);
        }
    }
    //pegando os valores da matriz B
    for(i=0;i<J_ORDER;i++){
        fscanf(pFile, "%f", &B[i]);
    }
	
    //pegando o resultado esperado para o teste do J_ROW_TEST
    resultadoEsperado = B[J_ROW_TEST];
	
    //primeira iteracao
	k=0;
    for(i=0;i<J_ORDER;i++){
        for(j=0;j<J_ORDER;j++){
            if(i!=j)
                novaMA[i][j]=MA[i][j]/MA[i][i];//Atualizando MA[i][j] com i!=j
            else
                novaMA[i][j]=0;//zerando a diagonal
			k=(i*J_ORDER)+j;
			novaMA2[k] = novaMA[i][j];
        }
        B[i]=B[i]/MA[i][i];//Atualizando B[i]
        X[i]=B[i];//pegando o primeiro valor de X[i]
    }
	
	//alocando memoria na GPU
	cudaMalloc((void**)&dev_novaMA, J_ORDER*J_ORDER*sizeof(float));
	cudaMalloc((void**)&dev_X, J_ORDER*sizeof(float));
	cudaMalloc((void**)&dev_novoX, J_ORDER*sizeof(float));
	cudaMalloc((void**)&dev_B, J_ORDER*sizeof(float));
	cudaMalloc((void**)&dev_maxDif, J_ORDER*sizeof(float));
	cudaMalloc((void**)&dev_maxnovoX, J_ORDER*sizeof(float));
	
	
	//copia os vetores X, B, novoX e novaMA para a GPU
	cudaMemcpy(dev_novaMA, novaMA2, J_ORDER*J_ORDER*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(dev_B, B, J_ORDER*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(dev_X, X, J_ORDER*sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy(dev_novoX, novoX, J_ORDER*sizeof(float), cudaMemcpyHostToDevice);
	
	/*pegar o tempo inicial*/
    //time(&time_start);
	clock_t start = clock(); 
	

    //testa convergencia do sistema e parte para as iteracoes
    if(convergenciaLinha(novaMA, J_ORDER) || convergenciaColuna(novaMA, J_ORDER)){//testes de convergencia
		
		//calcx(J_ORDER, dev_novaMA, dev_B, dev_X, dev_novoX, dev_maxDif, dev_maxnovoX);

		resolve<<<J_ORDER,1>>>(J_ORDER, dev_novaMA, dev_B, dev_X, dev_novoX, J_ERROR, J_ITE_MAX, dev_maxDif, dev_maxnovoX);

		//atualiza o valor de X no host
		cudaMemcpy(X, dev_X, J_ORDER*sizeof(float), cudaMemcpyDeviceToHost);

        //teste da fila a ser avaliada J_ROW_TEST
        for(i=0;i<J_ORDER;i++){
            resultado = resultado + MA[J_ROW_TEST][i]*X[i];
        }

        //imprimindo resultados na tela, como esta no exemplo da descricao do trabalho
        //comparacao do resultado obtido com o resultado esperado
        printf("\n\nIterations: %i\nRowTest %i => [%f] =? %f\n\n", iteracoes, J_ROW_TEST, resultado, resultadoEsperado);
    }//fim do metodo, caso o sistema passe em algum teste de convergencia

    //caso nao passe em nenhum teste de convergencia
    else
        printf("\nNao passou em nenhum dos testes de convergencia... encerrando programa...");

    /*pegar o tempo de termino*/
    //time(&time_end);

    /*fazendo a diferenca para pegar o tempo corrido*/
    //time_dif = difftime(time_end, time_start);

    /*imprime tempo corrido*/
    //printf("Tempo de execução = %.4f segundos\n", time_dif);
	printf("Time elapsed: %f\n", ((double)clock() - start) / CLOCKS_PER_SEC); 

    //liberando o espaço alocado na CPU
    free(MA);
    free(novaMA);
    free(B);
    free(X);
    free(novoX);
	free(novaMA2);
	free(vmaxDif);
	free(vmaxnovoX);
	
	//liberando espaco alocado na GPU
	cudaFree(dev_novaMA);
	cudaFree(dev_B);
	cudaFree(dev_X);
	cudaFree(dev_novoX);
	cudaFree(dev_maxDif);
	cudaFree(dev_maxnovoX);

    return 0;
}

//calcula o valor absoluto de "valor"
float modulo(float valor){
    if(valor<0)
        valor = valor * (-1.0);
return valor;
}

//funcao para verificar se o sistema converge pelo criterio das linhas
int convergenciaLinha(float **novoMA, int J_ORDER){
    int i,j;
    float soma;
    float maximo=0;
    for(i=0;i<J_ORDER;i++){
        soma = 0;
        for(j=0;j<J_ORDER;j++){//soma os valores absolutos da linha i
            soma = soma + modulo(novoMA[i][j]);
        }
        if(soma>maximo)//atualiza a maior soma
            maximo = soma;
    }
    if(maximo < 1)
        return 1;
    else
        return 0;
}

//funcao para verificar se o sistema converge pelo criterio das colunas
int convergenciaColuna(float **novoMA, int J_ORDER){
    int i,j;
    float soma;
    float maximo=0;
    for(j=0;j<J_ORDER;j++){
        soma = 0;
        for(i=0;i<J_ORDER;i++){//soma os valores absolutos da coluna j
            soma = soma + modulo(novoMA[i][j]);
        }
        if(soma>maximo)//atualiza a maior soma
            maximo = soma;
    }
    if(maximo < 1)
        return 1;
    else
        return 0;
}
